arXiv:1505.05524v3 [cond-mat.stat-mech] 13 Aug 2015 


Full self-similar solutions of the subsonic radiative heat equations 

Tomer Shussman^’^’|^ and Shay I. Heizler^’^’j^ 

^Raymond and Beverly Sackler School of Physics & Astronomy, 

Tel Aviv University, Tel Aviv 69978, ISRAEL 
^Department of Plasma Physics, Soreq Nuclear Research Center, Yavne 81800, ISRAEL 
^Department of Physics, Bar-Ilan University, Ramat-Gan, IL52900 ISRAEL 
‘^Department of Physics, Nuclear Research Center-Negev, 

P.O. Box 9001, Beer Sheva 84190, ISRAEL 

Abstract 

We study the phenomenon of diffusive radiative heat waves (Marshak waves) under general 
boundary conditions. In particular, we derive full analytic solutions for the subsonic case, that 
include both the ablation and the shock wave regions. Previous works in this regime, based on the 
work of [R. Pakula and R. Sigel, Phys. Fluids. 443, 28, 232 (1985)], present self-similar solutions 
for the ablation region alone, since in general, the shock region and the ablation region are not self¬ 
similar together. Analytic results for both regions were obtained only for the specific case in which 
the ratio between the ablation front velocity and the shock velocity is constant. In this work, we 
derive a full analytic solution for the whole problem in general boundary conditions. Our solution 
is composed of two different self-similar solutions, one for each region, that are patched at the heat 
front. The ablative region of the heat wave is solved in a manner similar to previous works. Then, 
the pressure at the front, which is derived from the ablative region solution, is taken as a boundary 
condition to the shock region, while the other boundary is described by Hugoniot relations. The 
solution is compared to full numerical simulations in several representative cases. The numerical 
and analytic results are found to agree within 1% in the ablation region, and within 2 — 5% in the 
shock region. This model allows better prediction of the physical behavior of radiation induced 
shock waves, and can be applied for high energy density physics experiments. 
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I. INTRODUCTION 


Radiation heat waves play important roles in many high energy density physics (HEDP) 
phenomena. In particular, they have major importance in inertial conhnement fusion (ICF) 
and in astrophysical and laboratory plasma [IHO]- In these experiments, laser beams deliver 
energy to the interior of a high-Z hohlraum that is converted into x-rays. Re-emission and 
further absorption of the x-rays in the cavity walls helps achieving a thermal source which 
acts as the drive for the experiments. The radiation is absorbed and contained within the 
cavity in a form of a radiative heat wave propagating through the hohlraum walls. It is 
therefore important to understand this phenomena, as a key to interpreting the experiments 
and the numerical simulations. 

The mechanism of the radiative heat waves is as follows: We consider a semi-inhnite wall, 
whose boundary is held at a high temperature. Usually, the temperature and density in this 
regime are such that the radiation energy is negligible compared to the matter energy, but 
the radiation heat flux is the dominant energy transport mechanism Dm The hot boundary 
radiates and heats the rest of the wall via photon transport. In the optically thick limit, 
a diffusive heat wave, characterized by a sharp temperature rise, propagates through the 
wall. If the wave propagates much faster than the speed of sound, hydrodynamic motion 
is negligible in the problem, and the wave is considered to be “supersonic” (Fig. 08)). 
If, however, the wave propagates slower than the speed of sound, the high matter pressure 
causes ablation of matter in the opposing direction of the heat wave. In addition, the 
heat wave is overtaken by a shock wave (Fig. 0b)). generated by the ablation pressure 
(from momentum conservation). The nature of the heat wave is temperature and density 
dependent, and can vary with time, as a supersonic diffusive front decelerates and becomes 
subsonic if the boundary temperature doesn’t change by much. If the temperature rises fast 
enough, the diffusive front accelerates and becomes more and more supersonic. In this work, 
we assume that the radiation and matter are in local thermodynamic equilibrium (LTE), 
which means they are strongly coupled. 

The hrst description of heat waves was proposed by Marshak [8], who obtained exact 
solutions for the radiative flow in the supersonic regime, in which hydrodynamic motion 
is negligible. For the subsonic case, in which hydrodynamics cannot be neglected, a full 
self-similar solution combining the ablation region and the shock region cannot be proposed 


2 


Subsonic 


Supersonic 




FIG. 1. (Color online) (a) A schematic diagram of a supersonic radiative heat wave. The heat wave 
propagates faster than the speed of sound and thus, hydrodynamic motion is negligible and the 
density is unchanged. The boundary condition on the surface is of a time dependent temperature. 
m is the Lagrangian coordinate, (b) A schematic diagram of a subsonic wave. Two separate regions 
exist: the ablation region, in which the density is low and the heat flux is dominant (left), and the 
shock region, in which the density is high and the heat flux is negligible (right). 

for the general case, since the problem is not self-similar altogether. Pakula and Sigel [MI], 
obtained self similar solutions of the heat wave in the subsonic case of an inhnitely dense wall, 
thus solving the ablation region only. Hammer and Rosen nails] proposed new solutions 
in this region, based on a perturbation theory. Self-similar solutions of the ablative subsonic 
regime in specihc cases were obtained in many other works [TUfTU] . In particular. Gamier et. 
al. [H] proposed a self-similar solution that includes both the ablative and shock regions, for 
a specihc case that ensures constant density over time (which is in this particular case, self¬ 
similar). These solutions are widely used for obtaining a better understanding of the heat 
wave phenomenon, evaluating the achieved temperature in IGF experiments [SHS] |20], or 
modeling hydrodynamical instabilities via linear perturbation amplihcation technique |21j . 

However, none of the previous works provides a full treatment of the shock wave and 
the ablation-shock interface for the general case. Only naive approximations, considering a 
constant ablation pressure were used to describe the shock region [I91I22H25]. These approx¬ 
imations are inaccurate in many cases, where the ablation pressure varies signihcantly over 
time. In this work, we propose a complete solution for both parts of the heat wave. We solve 
each part separately, in a self-similar fashion, and hud a continuous way to mathematically 
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patch them. For the ablative region we follow the solution of Ref. [9] and expand it, hnding 
a relation between the units of the conserved physical quantity and the temporal behavior of 
the boundary condition. Furthermore, we derive explicit expressions for the heat wave front 
coordinate and the energy contained within the heat wave, for general boundary conditions. 
In the shock region, we obtain a self-similar solution assuming a time-dependent pressure 
boundary condition. Then, the pressure at the ablation front, as obtained in the ablation 
region solution, is used as the boundary condition of the shock solution. The full solution 
is composed of these self-similar solutions. The results agree with full numerical simulation 
of the problem. 

This paper proceeds as follows: In Secs. [IT] and we solve the supersonic and the 


subsonic heat wave in the ablation region, using the mechanism of [9]. In Sec. IV we present 
the shock wave solution. In Sec. |V] we patch the solutions, provide numerical simulation 
results, and compare them with the integrated solution. A short discussion is presented in 
Sec. m 


II. SUPERSONIC WAVES 


A. Statement of the problem 


We consider a semi-inhnite wall of matter at density po- At time to fhe wall-vacuum 
interface is brought into contact with a thermal bath whose temperature is T{t). If hy¬ 
drodynamic motion is negligible, the radiative heat transport is described by one equation 
alone [DEIES]: 


de Of cIr d . ^ 


( 1 ) 


Where Ir = ^/krp is the Rossland mean free path, e is the internal energy per unit mass, 
and a = is the radiation density constant {a is Stefan-Boltzmann constant) and c is the 
speed of light. We assume the opacity and internal energy of the matter can be expressed 
in the form of power laws of the density and temperature, and follow the notation of 


— = srV 

I^R 


(2a) 


e = /ry-a 


(2b) 
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Using this assumption, Eq. becomes 


" " dm) 


(3) 


dt dm 

when A = 3f/3pQ'^~^^/16ag is a dimensional constant typical to this problem, and m = PqX 
is the Lagrangian coordinate. If we further assume that the boundary temperature is given 
as a power law, 

T{t) = Tot\ (4) 


then the problem has three typical units, mass [M], length [L] and time [6], and is charac¬ 
terized by exactly three parameters: 

[t] = [6] (5a) 

[T,] = [L]l[e]-h^ (5b) 

o 8+20+2/3 -2+3S-8 

[A] = [M]-^[L]-^[e]-i^ (5c) 

This means that a self-similarity of the hrst kind exists [T], and a construction of any 
dimensional variable using these parameters is unique. In this analysis, we distinguish 
between the two variables (m, t) of the problem and its three physical units. Although the 
problem is two-dimensional in the usual manner, it possesses three different physical units. 
We note that the temperature is related to the internal energy through Eq. [2bl so the units 
of the temperature can be dehned by forcing fP q^ to be dimensionless. The problem can 
also be solved using four units and / as an additional dimensional constant. 


B. The self-similar equation 


Using the parameters and the dimensions given above, a dimensionless parameter connect¬ 
ing the Lagrangian coordinate (with dimensions of [M]^[L]“^) and the temporal coordinate 
can be constructed: 


1 ^ a. A x(^ —a — 4 ) — 1 

0 


^ = mA^T^ ^ t 2 


The temperature prohle through the wall is therefore given as: 


( 6 ) 


Eq. [^also implies 


T(m,f)=ToUT(e) 

df _ i df 

dm m d^ 


(7) 


(8a) 
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dT 


v^dT 

t 


(8b) 


Here, = [t(/? — a — 4) — 1]/2 is the temporal power of Substituting these relations into 
Eq. [^yields a dimensionless Ordinary differential equation (ODE) 




—a—4 


dT 


dT 


+ = (« + 3)T-' hfE + 






d'^T 

W 


(9) 


The boundary conditions are T{^f) = 0, —)■ —oo and the parameter is determined 

uniquely by the normalization condition T(0) = 1. The total energy is given by Eq. 


E(t) = e(p, T)dm = f' 

Jo S-F Jo 

4+Q+/3 Tf4+Q+5') 1 ^ 

fP^A-2T„ » i^^+i / Tl>(i)d^ 

Jo 

The albedo can be determined from the relation between the absorbed energy and the 
emitted flux 


1 — a E{t) 

a aTot^'^ 


( 11 ) 


C. Boundary condition and conserved quantities 


For a given material, with known a, jj, the value of t fully determines the self-similar 
solution of Eq. The parameter t is strongly related to the units of the conserved quantity 
of the problem. If as a conserved quantity of dimensions 

[K] = [M]^^[L]^^[e]^^, ( 12 ) 


the self-similarity of the problem assures us that a dimensionless constant relates K to the 
other dimensional parameters. A, Tq, and t. The conservation of K means that this relation 
is independent of t: 

i = KA^^T^^ (13) 


Using Eq. [^we deduce: 


Wi = 


Ai 


W2 = 


/ 3 X 2 cr -|- /3 -|- 4 


Ai 


(14a) 

(14b) 
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2 ^(2a + 8-3,3)-2A3 

T = - — - ( 14cl 

j3 ( 3 X 2 + (a + /? + 4)Ai 

Similarly to the dimensionless constant ^ is determined by the normalization of the 
temperatnre prohle and K. For example, the case of constant net heat flnx throngh the 
bonndary, Sq is characterized by dimensions 

[So] = [MY[Lne]-^ (15) 

This yields t = and 

_ 1 i+a+j3 

So = ^ (16) 

The physical meaning of Sq is that the total energy obeys E{t) = Sot. Substitnting this in 
Eq. [ 1 ^ yields the constant value of ^ 

I =/p-" r ?'’(?)<« ( 17 ) 

Jo 


D. Solution of the equation 

For solving the equations, one must hnd for which T(0) = 1. This can be done 
using a shooting method, or using the self-similar coordinate relation (see Appendix 0. 
Once solved, the numeric value of ^p and the self-similar prohle can be used to obtain 
quantitative expressions for the heat front Lagrangian coordinate and the total absorbed 
energy per unit area, for given surface temperature and time. The expressions are of the 
form: 


M-A i-0 + a T(4-/3+a) + l 

mp = mop '2 Tq ^ t 2 


(18a) 


E 


_(£+A i+a+13 

eoPo ^ T, - t^- 


1 

2 


( 18 b) 


As an example, we take a medium of An, and use the values shown in na for the opacity 
and Equation of state (EOS) of the material. The values are specihed in Table Solving 
for the case of constant surface temperature yields the results: 


mp = 11.53 ■ IQ-Vo 


' g ' 
-cm^- 


E 


0.29po°-^^ 


2^3.55^0.5 


hJ 

mm2 


(19a) 

(19b) 
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Where Tq is measured in HeV and t is measured in nsec. For the case of constant boundary 
absorbed heat flux, which is important for hohlraum energy balance analysis [31 123], the 
temperature obeys T{t) = and we obtain the results: 


-4„-0.03Tnl.95j.0.775 


mF = 8.79 ■ 10 Po 

E = 0.21po°-^^To^-®® 




g 

cm^ 


hJ 


mm^ 


(20a) 

(20b) 


Quantitative expressions for the general boundary conditions are given in Fig. Numer¬ 
ical simulations (which will be presented in Sec. IV) yield the same expressions to within 
1% accuracy. 


TABLE I. Power law fits for the opacity and EOS of An in temperatnres 1 — 3HeV [ 12] 


Physical Quantity 

Numerical Value 

/ 

3.4 [MJ/g] 

(3 

1.6 

h 

0.14 

g 

1/7200 [g/cm2] 

a 

1.5 

A 

0.2 

r = (7 - 1 ) 

0.25 


III. SUBSONIC WAVES 

A. Statement of the problem 

In the subsonic case, the speed of sound exceeds the heat front velocity, and the hydro- 
dynamic motion is not negligible. The ablated matter density is much lower than the initial 
bulk density, while at the heat wave front the matter density is high enough to halt the 
thermal heat conduction. At this point, the high ablation pressure drives a shock through 
the bulk. In order to calculate the behavior of the system, one must solve the full radiative 
hydrodynamic equations (conservation of mass, momentum, and energy): 

dV du 
dt dm 

8 


0 


(21a) 



















FIG. 2. (Color online) (a) the constant of the Lagrangian coordinate given by Eq. [l^a). In circles, 
are marked the two special cases of constant boundary temperature (Eq. [I^a)) and constant 
absorbed heat flux (Eq. [20)(a)). (b) The constant of the energy per unit surface given by Eq. 
18 b). In circles, are marked the two special cases of constant boundary temperature (Eq. [I^b)) 


and constant absorbed heat flux (Eq. 20 b)). 


du 

m 


dP 

dm 

d 


= 0 

c d{aT^) 


de dV 

-h P — = 

dt dt dm dm 


(21b) 

(21c) 


Here, V = Yp is the specific volume, u is the matter velocity, P is the pressure, and m{x, t) = 
^ p{x', t)dx' is the Lagrangian coordinate. We assume the heat capacity and Rossland mean 
opacity follow Eq. and the EOS is well described by an ideal gas: 


P(p,T) = rpe(p,T) = - l)pe(p,T) 


( 22 ) 


From these relations we solve for the temperature: 

T = ‘ 


rf 


Substituting Eqs. ^ and ^ in Eq. |21c yields: 

IdPV^ + ^ (v^— (PV^-^^)^ 

r dt dt dm \ dm ^ ' 

Here, P is a dimensional parameter which is dehned as: 

16cr 


B = 


3(4 + 


, -* + °‘ 

3(Tf) t 


(23) 


(24) 


(25) 
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The problem again has three dimensions, and is characterized by three dimensional param¬ 
eters: 


W = [^] 

[To] = [T]^[M]t[0]-|-^ 

2^+A/3-m{4+q:) 3/j(4+a)-2Q;-2/3-3A/3-8 8-3/3+2a 

[B] = [M] - 3 - [L] - 3 - [0]^^ 


(26a) 

(26b) 

(26c) 


In this case we dehne / to be dimensionless (instead of /po ^ in the snpersonic case), and thus 


the units of Tq are now different from Eq. We note that there exists another dimensional 
parameter in the problem, the initial density po- This parameter can be neglected in the 
highly subsonic regime, and in the ablation region only, due to the fact that the matter 
density in the ablation region is much lower than po. Therefore, it can be assumed that the 
density goes from 0 at the rear surface and approaches inhnity in the ablation front mm- 
Neglecting the initial density of the problem, prevents us from solving the shock region 
and the ablation region altogether. In fact, the full problem, which includes both the 
ablation and shock is not self-similar, since the shock region depends on the initial density 
of the matter, via Hugoniot relations, while the ablation region depends on the dimensional 


parameter B due to the heat flux. Therefore, solving Eqs. while neglecting the initial 
density yields a solution for the ablation region alone, in which the temperature at the 
ablation front is zero (instead of the shock temperature), and the front density approaches 
inhnity. The same discussion can be applied for the how velocity. For solving the ablation 
region, we must neglect the how velocity at the heat front, and assume it approaches 0 
We will now solve the ablation region, following a method similar to the one used in [HI 

[nmn. 


B. The self-similar equations 

Every dimensional variable can be parameterized as a power law of the dimensional 
parameters: 

X = (27) 

In addition, the self-similar coordinate is parameterized as following: 

^ ( 28 ) 
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The powers wxi and Wi are deduced from equations Specifically, the Lagrangian coordi¬ 
nate, with units of [M]/[L]^ is given: 




1 

^ ^j^At-22i2/3-2a-8-/3A+(4+a)At^-2-2(4+Q:-/3)T+M(3+(4+a)T)-A(2+/3T)j 4+2A-4 m 


(29) 


The power law dependence of the rest of the dimensional variables is specihed in Table [ITj 


TABLE II. Power law dependence of the physical quantities in the subsonic ablative 
heat solution 


Physical variable 

Power-law dependency 

Lagrangian 

Coordinate 

^ ^ ^^2—2/3-I-/3A—(4+o)/4 

^2-|-2(4+a-/3T-/2(3+(4+a)T)-|-A(2-|-^T)l 4+2A-4 m 

Temperature 

II 

rs- 

Velocity 

~ /^-/42V(2+A)-(4+o)/i^^_|_^(2 + A)T-(4-|-Q:)/4T'j 4+2A-4M 

Pressure 

^ ^ 1 I 

PljTl t) _ P 2-\-X—2fj. 

Specific Volume 

V{m,t) = V |^5-l2^-4-«+2/3^1-(4+«-2/3)T^ 2+^-2/- 

Heat Flux 

S{m, t) = S (^52-3/4238+2a+2/3+3A/3-3(4+a)/4 

^-2+3/i+(2(4+a+/3)+3,9A)T-3(4+a)/4T^ 4+2A-4 m 

Energy 

E{m, t)=E (^ 52 - 3 M 7 ^ 8 + 2 «+ 2 / 3 + 3 A/ 9 - 3 ( 4 +a)M 

^2+2A-At-|-(2(4+a+/3)-|-3/3A)T-3(4+a)/2T^ 4+2A-4 m 


Using Table and the equivalent of Relations for the subsonic case, a set of dimen¬ 
sionless ODEs can be obtained: 


d \ ~ du 

U.V3 + «.3f5^jr--=0 


(30a) 


. 9 \ _ dP 

u-H — = 0 
dij di 


(30b) 
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dP ~dV 

{WV3 + WP3)PV + W3^ ( + P-^ 


+ P ( Wv3 + W3^-^ ) V = 




(9 1 d( 




+ 


A I O 4+a —2/3 

4 + « - /^T~A P 






+ 


yA j 


4+a —/3 




n-a)PvP^ + v'-'‘^^ 


(30c) 


The boundary conditions are = 0, = 0, +(0) = 0, ^|^=o = 0, m(^f) = 0, and 

the free parameter is determined from the normalization condition of the temperature, in 
a manner similar to the supersonic solution. We note that since the temperature appears in 
Eq. [^only via the relation specified in Eq. ^ the normalized solution should now satisfy 
T(0)^ = P(0)E(0)^“^ = 1. Then, the absorbed energy is given by: 


./n T 2 ./n T 2 


C. Boundary condition and the conserved quantities 

The self-similarity of the problem assures us that a dimensionless constant relates K to 
the other dimensional parameters, B, Tq, and t. The conservation of K means that this 
relation is independent of t: 

i = (32) 


Using Eq. 26 we deduce: 


[2/3 P A/3 — /i(4 -|- a)] wi -I- fiW 2 — /3Ai 

[3/i(4 + a) — 2a — 2(3 — 3A/3 — 8] tci -|- (2 — 3 fi)w2 = (3X2 

[rci(8 — 3/3 -f 2a) -|- / 3 A 3 ] 2 

I3w2 (3 


(33a) 

(33b) 

(33c) 


12 
























Specifically, for the case of constant flnx, with nnits [S'o] = [M][6*] the bonndary 
condition obeys: 


2 - 3/i 

8 -|- 2q; -p 2/3 + 3/3A — 3(4 -1- cx)^ 


(34) 


D. Solution of the equations 


The numerical solution of Eqs. is obtained by a double shooting method, or by using 
the self-similar relation (see Appendix]^. The self-similar prohles for the case of t = 0, as 
obtained by the numerical integration of Eqs. ^ are presented in Fig. Once solving the 


dimensionless quantities, we can insert them to the self-similar relations given in Table [IT] 
and obtain quantitative expressions for the heat front Lagrangian coordinate, total energy, 
and the ablation pressure at the front. For the case of constant surface temperature: 


mp = 10.17- 

[mm"* 

Pf = [Mbar] 


Lcm^ j 


For the case of constant boundary net-flux we obtain T{t) = and: 


(35a) 

(35b) 

(35c) 


mp = 8.32 ■ 


Lcm^ 


E = 0A2l-^h 


hJ 


mm^ 


Pp = 


2.63+-0.125 


[Mbar] 


(36a) 

(36b) 

(36c) 


Quantitative expressions for the general boundary conditions are given in Fig. Here, 
again, numerical simulations reproduce the analytic expressions to within 1% accuracy. 
In addition, it reproduced both the analytic and numeric solutions of the supersonic and 
subsonic ablative heat waves, which were obtained and checked in previous works MEol 
[23]. 
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FIG. 3. (Color online) The dimensionless variables, obtained in the numerical integration of Eqs. 


30 for the case of constant temperature (t = 0). Since the shock region is neglected, the temper¬ 


ature drops sharply to 0 at and the density diverges. The ablation pressure (e.g. the pressure 
at the ablation front) has a finite value. 


IV. SHOCK SOLUTION 


In this section we solve the shock region of the subsonic heat wave. As mentioned earlier, 
the full problem, composed of the shock region and the ablation region, is not self-similar. 
Therefore, we solve the shock region independently, and in the next section we patch both 
regions and obtain a full solution. 


A. Statement of the problem 

Mathematically, the shock region is different from the ablation region by two main pa¬ 
rameters. Since the downstream and upstream densities are related via Hugoniot Relations, 
the initial density po is a key parameter of the shock region, and it seems as if the shock 
region is not self-similar at all. However, in the case of a strong shock, as obtained in highly- 
subsonic heat waves, radiation heat conduction is negligible in the shock region. Therefore 
we can ignore the dimensional parameter B, which is presented in the previous section. 

The question remains what is the proper boundary condition to describe the surface 
between the regions, and three natural contenders are the ablation pressure, the ablation 
velocity, and the front density. The density is obviously a bad choice since it diverges at the 
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FIG. 4. (Color online) (a) the constant of the Lagrangian coordinate given by Eq. 28 In circles, 


are marked the two special cases of constant boundary temperature (Eq. [3^a)) and constant 


absorbed heat flux (Eq. ^ma)). (b) the constant of the energy per unit surface given by Eq. 27 


In circles, are marked the two special cases of constant boundary temperature (Eq. 35 b)) and 


constant absorbed heat flux (Eq. 36 b)). (c) The same for the ablation pressure. 


ablation front, and the velocity is assumed to be zero at the front. In fact, matching the 
self-similar solution in the ablation region to a self-similar solution in the shock region will 
only be valid if the flow velocities in the shocked region are much smaller than the exhaust 
flow in the ablation region, which is of order the speed of sound of the radiatively-heated 
material. 

Regarding the discussion above, the pressure boundary condition should be satisfying. 
We use the heat front pressure as a boundary condition, and assume that the shock boundary 
is well approximated as a surface with power law pressure dependence. 
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(37a) 


For the shock region in the general case, the set of equations needs to be solved is: 

dV du ^ 
dt dm 
du dP 


dt dm 


We assume the boundary condition: 


P{t) = Pof^ 


(37b) 

(37c) 

(38) 


In the case of a shock created by a heat wave, this boundary condition is obtained from 
Table [TT} We further assume that the shocked material EOS obeys the ideal gas relation (Eq. 
22) where r is not necessarily the same for the ablation region and the shock region. In the 
shock front, the boundary conditions obey Hugoniot Relations [1]: 


D — u D 

Vs ^ Vo 


Ps = 



PsUsVo = D{es+ " 1 ^ 


(39a) 

(39b) 

(39c) 


Here, Vq is the unperturbed specihc volume. Vs, us, Psi and es are the specihc volume, 
matter velocity, pressure, and thermal energy right after the shock, and D is the shock 
velocity in the lab frame. The parameters of the problem and their dimensions are: 


W = [^] (40a) 

[Po] = [L]-^[M][9]-^-^^ (40b) 

[W] = [Lf[M]-^ (40c) 

Using these parameters, we again deduce self-similar relations of the form: 

X = XPo""^ (41) 


The power law dependence of the physical variables is given 
self-similar coordinate is: 


in Table III 


Specihcally, the 


^ = m(^ 




(42) 
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TABLE III. Power law dependence of the physical quantities in the shock solution 


Physical variable 

Power-law dependency 

Lagrangian Coordinate 


Velocity 

u{m, t) = uPq Vq 

Pressure 

P{m, t) = PPqP^ 

Specific Volume 

II 

Energy 

E{m,t) = EP^vjp+l^s 


B. The Self-Similar equations 


Substituting Eq. 
. ^ A 

dm m 


42 


into Eq. 37, and using the derivatives % = — (1 +'^s/ 2 ) | and 


^ , we obtain the self-similar set of equations: 

rim. m ' tl 


+ -- = 0 

xs\ du dP 


2 J r di r 




(43a) 

(43b) 

(43c) 


The boundary conditions are T’(O) = 1 and the self-similar Hugoniot Relations at ^s- 
The shock velocity obeys: 


D = Vom = (^ + y) ^ 

Substituting Eq. M in Eq. yields the self-similar Hugoniot relations: 

Us 


Ks = l- 


(! + ¥)« 


= (1 + y) 


1 + y ) & 


Rearranging the equations yields: 

y = 




7-1 


r + 2 


7 + 1 


(44) 

(45a) 

(45b) 

(45c) 

(46a) 
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(46b) 


s* = ^ (l + y') 

p _ul r 

-TC — —— 


2 


n ^ + 2~ 
D = -^^us 


(46c) 


(46d) 


Once again, the free parameter is determined from the normalization condition at 
^ = 0, bnt this time the bonndary pressnre is normalized to match the bonndary condition. 

C. Solution of the Equations 


The set of Eqs. 43 is solved using a shooting method, just like in the supersonic case. After 


obtaining the self-similar prohles, quantitative expressions for the shock front Lagrangian 


Coordinate and the bulk velocity are derived, using Table III (The shock velocity can be 


calculated using Eq. 46d). For convenience, we choose the same value of r that we used for 


the ablation region. For a constant boundary pressure we get: 


ms = 4.66 ■ 


cm^ 


Us = 2.15P; 


0.5 


km 


sec 


(47a) 

(47b) 


The cases of constant surface temperature in the ablation region, and constant absorbed 
flux are of particular interest for this work. In the case of constant temperature, the ablation 


pressure temporal behavior yields (see Eq. 35) wp^ = Xs = —0.45 (we remind that the 

as a 


notation of wp- is dehned in Eq. 27 and calculated in Table [^. Substituting this 
boundary condition for the shock gives the relation: 


-3 dO.5+0.7765 


ms = 7.34 ■ lO-^Pn f 


30.5+-0.2235 


cm^ 
km 


sec 


(48a) 


(48b) 


Us = 2 . 62 P^-H- 

In the case of constant absorbed flux, the ablation pressure behaves as wp^ = Xs = 
-0.124, and the numerical expressions are: 

g 


ms = 5.17-10 


-3 p0.5j.0.938 
-in 




cm^ 


Us = 2.23P^-H 


0.5+-0.062 


km 


sec 


(49a) 

(49b) 
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FIG. 5. (Color online) (a) the constant of the Lagrangian coordinate given by Eq. 42 In circles, are 
marked three special cases of constant boundary temperature (Eq. [4^a) ), constant absorbed heat 
flux (Eq. [4^a)) and constant ablation pressure a)), (b) the constant of the matter velocity at 
the shock front. In circles, are marked three special cases of constant boundary temperature (Eq. 


48 b)), constant absorbed heat flux (Eq. 49 b)) and constant ablation pressure (47 b)). 


Quantitative expressions for the general boundary conditions are given in Fig. Using 
these expressions, one can derive the shock solution for any given ablation region boundary 
condition. This can be done by substituting Tg with wp^ and Pq, shock with Po,heat7)7^% where 
Po,shock, Po,heat correspoud to the shock region boundary condition constant, and the ablation 
pressure constant, respectively. We note that the temperature can be obtained, using Eq. 


23 In this work, for simplicity we assume the same values of /, /?, /i and r for both the 


ablation regime and the shock regime, although this is usually not the case. 


V. FULL SOLUTION AND NUMERICAL SIMULATIONS 


In this final section, we present a full solution of the subsonic heat equation, decomposed 


of the two self-similar solutions obtained earlier for the ablation region (Sec. Ill) and for the 


shock region (Sec. IV). The two solutions are patched together in the following manner: 


• The ablation region is solved using a specihed boundary condition, given in the form 
of Eq. 1^ Self-similar and quantitative profiles for the temperature, pressure, density 
and velocity are obtained up to the ablation front mp. 
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FIG. 6. (Color online) An example of the patching method of the two different regions, in a 
constant surface temperature problem. The two self-similar solutions of the two different regions 
are marked in dashed line, while the full solution is taken as the ablation solution from m = 0 to 
the patching point (when the curves collide), and from there, the shock solution is used. The small 
box in the velocity graph is a zoomed look on the shock region. 


• The ablation pressure is applied as the power law boundary condition of the shock 
region problem (Pq, shock = Po,heatPo’^^^^^s = wps). 

• The shock region is solved using this boundary condition and Hugoniot relations at 
the front. Self-similar and quantitative profiles for the temperature, pressure, density 
and velocity are obtained from m = 0 to the shock front ms > mp. 

• The full solution is given by the ablation region profile from m = 0 to m = mp, the 
shock region profile from m = mp to m = mg and the undisturbed matter profile from 
there forward. 
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FIG. 7. (Color online) A comparison between the self-similar shock solution and numerical sim¬ 
ulations. The dashed (full) line shows the shock analytic (numerical) solution for a boundary 
condition of P = 2.71t“^'^^Mbar. The dashed-dotted line is yielded from a full heat wave numer¬ 
ical solution in the boundary condition T = IHeV (this heat wave boundary condition yields an 
ablation pressure of P = 2.71t“‘^'^^Mbar exactly). The shock numerical solution agrees with the 
analytic solution within 1%, and the full simulation agrees with the analytic solution to within 3% 
in the shock region. 

An example of the patching method is shown in Fig. for the case of constant temper¬ 
ature boundary condition Tq = IHeV, and at the time t = 0.05ns. The solution is shown at 
an early time so that the full prohle can be well seen, as the shock wave propagates ahead 
of the heat front at later times. 

In order to solve for the self-similar shock region prohle, we had to assume that the 
boundary condition is well described as a surface with power law pressure dependence. This 
is obviously not the case, since the ablation front matter velocity is actually not zero, and 
the density does not diverge, as opposed to the assumptions that led to the analytic ablation 
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FIG. 8. (Color online) A comparison between the model (dashed lines) and simnlations (full lines) 
under a boundary condition of constant surface temperature (t = 0) with Tq = IHeV in early 
times. Presented are the profiles of the temperature, density, pressure and matter velocity. The 
small box in the velocity graph is a zoomed look on the shock region. 

region solution. However, comparison of the self-similar solution with the numerical solution 
of a shock, in the boundary condition of Pq = 2.71 and Xs = —0.447 (which corresponds 
to a heat wave of IHeV constant temperature), and with a full numerical simulation of a 
subsonic heat wave, as shown in Fig. shows an agreement of 1% between all three prohles. 
This is due to the fact that only the boundary of the problem is affected by the velocity 
condition, as the rest of the bulk is solely affected by the pressure gradient. 

For the numerical simulations we used a standard ID hydrodynamic Lagrangian code, 
which includes radiative transport in ID LTE diffusion approximation. The code was val¬ 
idated via regular test cases. In addition, it reproduced the analytic solutions of the su¬ 
personic and subsonic ablative heat waves, which were obtained and checked in previous 
works [MlEQlES]. In the simulations, we used the same power laws for the EOS and the 
opacity as in the self-similar model. In the “full simulation”, we set the surface bound- 
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FIG. 9. (Color online) A comparison between the model (dashed lines) and simulations (full lines) 
under a boundary condition of constant absorbed flux (Sq, t = 0.123) with Tq = IHeV in early 
times. Presented are the profiles of the temperature, density, pressure and matter velocity. The 
small box in the velocity graph is a zoomed look on the shock region. 


ary condition to be Tg{t) = for different values of t, and checked that the conserved 
quantities which correspond to each value of t were indeed conserved. 


In Figs. [8 IT the full analytic solution is compared to numerical simulation results in 
different boundary conditions and times. We present three representative cases: constant 
temperature, constant absorbed flux which yields t = 0.123, and constant ablation pressure 
which yields t = 0.17. In Figs. [8 W we present the comparison in relatively early times 
(0.05 ^ t ^ 0.15nsec). There is a good match between the analytic solution and the full 
numerical simulations, as the results agree to within 1% in the ablation region, and in the 
shock region, the agreement is about 5%. 


The lack of perfect agreement between the full solution and the analytic one, may be due 
to the fact that in early times, the heat front is not subsonic (and the ablation density is 
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FIG. 10. (Color online) A comparison between the model (dashed lines) and simulations (full lines) 
under a boundary condition of constant ablation pressure (Pqj 'Ts = 0, t = 0.17) with Tq = IHeV 
in early times. Presented are the profiles of the temperature, density, pressure and matter velocity. 
The small box in the velocity graph is a zoomed look on the shock region. 


not infinite). The transition from supersonic to subsonic also invalidates the approximation 
of negligible flow velocities, until the subsonic flow becomes established. 


In Fig. H we present a comparison in a later time, t = Insec. The match between the 
model and the simulation is even better, as in the shock region they agree to within 2 — 3%. 
The better match at late times is due to the build-up time of the solution. We note that the 
numeric solution front is and not at mp + ms as could naively be expected. We assume 
that the shock front is weakly affected by the ablation front position, due to causal delay. 

In addition, we see that the physical properties obtained using the naive approximation 
of constant ablation pressure, signihcantly differ from the exact full solution (for example, 
see the density prohle and shock front in Fig. 0. This demonstrates the major benehts of 
using the model presented in this work, for evaluating the physical properties of the shock 
region. 
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FIG. 11. (Color online) A comparison between the model (dashed lines) and simulations (full 
lines), under the different boundary conditions used in Figs. [8 10 in late times. Presented are the 
profiles of the temperature, density, pressure and matter velocity. The small box in the velocity 
graph is a zoomed look on the shock region. 


VI. DISCUSSION 

In this paper we investigated the diffusive radiative heat equation. At first we generalized 
the self-similar solutions of both the supersonic and the subsonic heat equations, for a general 
boundary condition of the type T = Toh^, and for power law opacity and heat capacity. We 
used the technique shown in [9] (that solved the equations for several specific cases) with 
the notation of na (i.e. Super-transition-arrays based An opacity). The numerical results 
of tuq and Cq for An are given in Fig. [^for the supersonic case and in Fig. |^for the subsonic 
case, and can be used for general purposes, such as hohlraum temperature evaluation in 

IGF [3HSlEni[23]. 

Next, we focused on finding a self-similar solution of the shock region, under the assump¬ 
tion of a given pressure boundary condition of the type P = PqP^ on one boundary of the 
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system, and the strong shock Hngoniot relations on the other. The solntion was fonnd to 
be in agreement with nnmerical simnlation (Fig. [^. 

Finally, We set the resnlting Pq and T5 that the ablation heat wave solntion yielded, 
snbstitnte them as inpnt to the shock solntion, and patched both solntions at the interface. 
Using this method, we obtained a fnll analytic solntion to the non self-similar problem, 
composed of two different solntion, each one of them is self-similar by itself and is valid in a 
different region of the problem. This simple techniqne was tested with nnmerical simnlations 
and was fonnd to be accnrate to within 3%. 

As mentioned in the Introdnction, the shock wave can be evalnated via a naive approxi¬ 
mation of constant ablation pressnre [22] , which is of conrse valid only if the ablation pressnre 
is relatively constant with time. Shock waves experiments have been analyzed nsing this 
approximation [231121] . In fntnre work we plan to reprodnce the ablative shock analysis, 
and check the agreement between this new model and experimental resnlts. 

Appendix A: A shooting method for solving the self-similar equations 

Obtaining a solntion for the self-similar eqnations reqnires solving a set of eqnations with 
bonndary conditions at both ends of the space. This is generally done using a shooting 
method, or using the self-similarity of the equation. In the case of the supersonic equation, 
Eq. 1^ one must hnd for which T(0) = 1. In both methods, is initially guessed. 

= 1 is usually a good estimate. In the shooting method, each step the value of ^ updates 
according to whether T(0) is larger or smaller than one. The step size decreases with the 
number of iterations. A newton shooting method [2] can easily be applied for this problem. 
The self-similar method for Ending ^p is based on the fact that for a given time and mass, 
the coordinate ^ will only depend on the boundary temperature, according to Eq. [^ This 
conservation of coordinates implies 

, I3-OC-4 

- - ( - ) ^ 

If T(0) 7^ 1, To in Eq. [^must be multiplied by a normalization factor of T(0). Therefore 
the relation between the normalized ^p and the guessed is 

p = 
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(A2) 





The subsonic set of equations can be solved in the same manner, using a double shooting 


method. For a given value of initial values for guessed and the set of 


9 ? 


Eqs. is numerically integrated using a simple ODE solver. Then, the guess is updated 
according to whether P(0) is larger or smaller then 0. In the latter case, the integration 
diverges at some .^ > 0. This procedure repeats for every iterative value while the 
determination of ^p is done in a manner similar to the supersonic solution, in order to 
normalize the self-similar quantities. 
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